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I. INTRODUCTION 



The prediction of ionic currents in protein channels of biological membranes is one of the 
central problems of computational molecular biophysics. None of the existing continuum 
descriptions of ionic permeation captures the rich phenomenology of the patch clamp ex- 
periments It is therefore necessary to resort to particle simulations of the permeation 
process Q-pj. Computer simulations are necessarily limited to a relatively small number 
of mobile ions, due to computational difficulties. Thus a reasonable simulation can describe 
only a small portion of the experimental setup of a patch clamp experiment, the channel and 
its immediate surroundings. The inclusion in the simulation of a part of the bath and its 
connection to the surrounding bath are necessary, because the conditions at the boundaries 
of the channel are unknown, while they are measurable in the bath, away from the channel. 

Thus the trajectories of the particles are described individually for each particle inside the 
simulation volume, while outside the simulation volume they can be described only by their 
statistical properties. It follows that on one side of the interface between the simulation 
and the surrounding bath the description of the particles is discrete, while a continuum 
description has to be used on the other side. This poses the fundamental question of how 
to describe the particle trajectories at the interface, which is the subject of this paper. 

We address this problem for Brownian dynamics (BD) simulations, connected to a prac- 
tically infinite surrounding bath by an interface that serves as both a source of particles 
that enter the simulation and an absorbing boundary for particles that leave the simulation. 
The interface is expected to reproduce the physical conditions that actually exist on the 
boundary of the simulated volume. These physical conditions are not merely the average 
electrostatic potential and local concentrations at the boundary of the volume, but also their 
fluctuation in time. It is important to recover the correct fluctuation, because the stochastic 
dynamics of ions in solution are nonlinear, due to the coupling between the electrostatic field 
and the motion of the mobile charges, so that averaged boundary conditions do not repro- 
duce correctly averaged nonlinear response. In a system of noninteracting particles incorrect 
fluctuation on the boundary may still produce the correct response outside a boundary layer 
in the simulation region jsj]. 

The boundary fluctuation consists of arrival of new particles from the bath into the 
simulation and of the recirculation of particles in and out of the simulation. The random 
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motion of the mobile charges brings about the fluctuation in both the concentrations and 
the electrostatic field. Since the simulation is confined to the volume inside the interface, 
the new and the recirculated particles have to be fed into the simulation by a source that 
imitates the influx across the interface. The interface does not represent any physical device 
that feeds trajectories back into the simulation, but is rather an imaginary wall, which the 
physical trajectories of the diffusing particles cross and recross any number of times. The 
efflux of simulated trajectories through the interface is seen in the simulation, however, 
the influx of new trajectories, which is the unidirectional flux (UF) of diffusion, has to be 
calculated so as to reproduce the physical conditions, as mentioned above. Thus the UF is 
the source strength of the influx, and also the number of trajectories that cross the interface 
in one direction per unit time. 

The mathematical problem of the UF begins with the description of diffusion by the 
diffusion equation. The diffusion equation (DE) is often considered to be an approximation 
of the Fokker-Planck equation (FPE) in the Smoluchowski limit of large damping. Both 
equations can be written as the conservation law 

The flux density J in the diffusion equation is given by 

J(x, t) = — [eVp(x, t) - f(x)p(x, t)} , (2) 

7 

where 7 is the friction coefficient (or dynamical viscosity), e = , ks is Boltzmann's 

m 

constant, T is absolute temperature, and m is the mass of the diffusing particle. The 
external acceleration field is f{x) and p{x,t) is the density (or probability density) of the 
particles jjj] . The flux density in the FPE is given by where the net probability flux density 
vector has the components 

J x (x,v,t) = vp(x,v,t) 

(3) 

J v (x,v,t) = - {jv - f(x))p(x,v,t) - e>y\7vp(x, v,t). 

The density p(x, t) in the diffusion equation with (0) is the probability density of the 
trajectories of the Smoluchowski stochastic differential equation 

x = ^(x) + ^w, (4) 



where w(t) is a vector of independent standard Wiener processes (Brownian motions). 

The density p(x, v, t) is the probability density of the phase space trajectories of the 
Langevin equation 

x + 7£c = f(x) + y/2ejw. (5) 

In practically all conservation laws of the type J is a net flux density vector. It is 
often necessary to split it into two unidirectional components across a given surface, such 
that the net flux J is their difference. Such splitting is pretty obvious in the FPE, because 
the velocity v at each point x tells the two UFs apart. Thus, in one dimension, 

poo rO 

JLn(x,t) = / vp(x,v,t)dv, J RL (x,t) = - vp(x,v,t)dv 

JO J -oo 

(6) 

POO 

J net( x ^) = J Ln{x,t) - J RL {x,t) = / vp(x,v,t)dv. 



In contrast, the net flux J(x,t) in the DE cannot be split this way, because velocity 
is not a state variable. Actually, the trajectories of a diffusion process do not have well 
defined velocities, because they are nowhere differentiable with probability 1 [10]. These 
trajectories cross and recross every point x infinitely many times in any time interval [t, t + 
At], giving rise to infinite UFs. However, the net diffusion flux is finite, as indicated in 



eq.(j2J). This phenomenon was discussed in detail in 11 1, where a path integral description 
of diffusion was used to define the UF. The unidirectional diffusion flux, however, is finite 
at absorbing boundaries, where the UF equals the net flux. The UFs measured in diffusion 
across biological membranes by using radioactive tracer bl are in effect UFs at absorbing 
boundaries, because the tracer is a separate ionic species [13]. 

An apparent paradox arises in the Smoluchowski approximation of the FPE by the DE, 
namely, the UF of the DE is infinite for all 7, while the UF of the FPE remains finite, even 
in the limit 7 — > 00, in which the solution of the DE is an approximation of that of the FPE 
jl^j ]. The "paradox" is resolved by a new derivation of the FPE for LD from the Wiener 
path integral. This derivation is different than the derivation of the DE or the Smoluchowski 
equation from the Wiener integral (see, e.g. Qj-j^) by the method of M. Kac j^|. The 
new derivation shows that the path integral definition of UF in diffusion, as first introduced 
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ll| . is consistent with that of UF in the FPE. However, the definition of flux involves 
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the limit At — > 0, that is, a time scale shorter than 1/7, on which the solution of the DE is 
not a valid approximation to that of the FPE. 

This discrepancy between the Einstein and the Langevin descriptions of the random 
motion of diffusing particles was hinted at by both Einstein and Smoluchowski. Einstein 
remarked that his diffusion theory is based on the assumption that the diffusing particles 
are observed intermittently at short time intervals, but not too short, while Smoluchowski 
[l5l | showed that the variance of the displacement of Langevin trajectories is quadratic in t 
for times much shorter than the relaxation time I/7, but is linear in t for times much longer 
that I/7, which is the same as in Einstein's theory of diffusion 

The infinite unidirectional diffusion flux imposes serious limitations on BD simulations 
of diffusion in a finite volume embedded in a much larger bath. Such simulations are used, 
or example, in simulations of ion permeation in protein channels of biological membranes 
l|. If parts of the bathing solutions on both sides of the membrane are to be included in 
the simulation, the UFs of particles into the simulation have to be calculated. Simulations 
with BD would lead to increasing influxes as the time step is refined. 

The method of resolution of the said "paradox" is based on the definition of the UF of the 
Langevin dynamics (LD) in terms of the Wiener path integral, analogous to its definition 
for the BD in [llj]. The UF JLn{x,t) is the probability per unit time At of trajectories that 
are on the left of x at time t and are on the right of x at time t + At. We show that the UF 
of BD coincides with that of LD if the time unit At in the definition of the unidirectional 
diffusion flux is exactly 

At = I (7) 

We find the strength of the source that ensures that a given concentration is maintained on 
the average at the interface in a BD simulation. The strength of the left source Jlr is to 
leading order independent of the efflux and depends on the concentration Cl, the damping 
coefficient 7, the temperature e, and the time step At, as given in eq.(|27j). To leading order 
it is 



Jlr = J^r-C L + 0[-\. (8) 



TvyAt \7 / 

We also show that the coordinate of a newly injected particle has the probability dis- 
tribution of the residual of the normal distribution. Our simulation results show that no 
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spurious boundary layers are formed with this scheme, while they are formed if new particles 
are injected at the boundary. The simulations also show that if the injection rate is fixed, 
there is depletion of the population as the time step is refined, but there is no depletion if 
the rate is changed according to eq.(|8|). 

In Section m we derive the FPE for the LD (j5J) from the Wiener path integral. In Section 
IIHI we define the unidirectional probability flux for LD by the path integral and show that is 
indeed given by ©• In Section HVl we use the results of (l3| to calculate explicitly the UF in 
the Smoluchowski approximation to the solution of the FPE and to recover the flux ©. In 

n 

Section El we use the results of (Tj| to evaluate the UF of the BD trajectories (jlj) in a finite 
time unit At. In the limit At — > the UF diverges, but if it is chosen as in (J7J), the UFs of 
LD and BD coincide. In Section IV1I we describe the a BD simulation of diffusion between 
fixed concentrations and give results of simulations. Finally, Section IVHI is a summary and 
discussion of the results. 

II. DERIVATION OF THE FOKKER-PLANCK EQUATION FROM A PATH IN- 
TEGRAL 

The LD (0) of a diffusing particle can be written as the phase space system 

x = v, v = — 7f + f(x) + w. (9) 

This means that in time At the dynamics progresses according to 

x(t + At) = x(t)+v(t)At + o(At) (10) 



v(t + At) = v(t) + \-iv{t) + f(x(t))]At + ^f2^| Aw + o(At), (11) 

where Aw ~ A/"(0, At), that is, Aw is normally distributed with mean and variance At. 
This means that the probability density function evolves according to the propagator 

Prob{a;(t + At) = x, v(t + At) = v} = p(x, v, t + At) = o(At) + (12) 
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To understand (jl2j) . we note that given that the displacement and velocity of the trajectory 
at time t are x{t) = £ and v(t) = r\, respectively, then according to eq- lfTUj) . the displacement 
of the particle at time t + At is deterministic, independent of the value of Aw, and is 
x = £ + V^t + o(At). Thus the probability density function (pdf) of the displacement is 
5(x — £ — rjAt + o(At)). It follows that the displacement contributes to the joint propagator 
(112}) of x(t) and v (t) a multiplicative factor of the Dirac 5 function. Similarly, eq. ()llj) means 
that the conditional pdf of the velocity at time t + At, given x(t) = £ and v(t) = r], is 
normal with mean i] + [—777 + /(£)]At + o(At) and variance 2e7 At + o(At), as reflected in 
the exponential factor of the propagator. If trajectories are terminated at the ends of an 
finite or infinite interval (a,b), the integration over the displacement variable extends only 
to that interval. 

The basis for our analysis of the UF is the following new derivation of the Fokker-Planck 
equation from eq. (|T2j) . Integration with respect to £ gives 

p(x,v,t + At) = o(At) + (13) 
1 r^-^t^lj'-'-^+fj'-MM' } dr). 



y/AejnAt J ^ j Ae^At 

Changing variables to 

v — 7] — [—7^ + f(x — r]At)}At 
~ U ~ y/2e^At ' 

and expanding in powers of At, the integral takes the form 

1 f°° 

p(x,v,t + At) = -= / e~ u2/2 dux (14) 



^7r(l -7At + o(At)) 



p(x - v(l + 7 At) At + o(At), v(l + 7At) + u^le^At - f{x)At{\ + 7 At) + o(At),t) 
Reexpanding in powers of At, we get 

p(x - v(l + jAt)At + o(At), v(l + jAt) + uy/2ei At - j{x)At{\ + 7 At) + o(At),t) = 

p(x, v, t) - vAt dp{ - X : V ^ + ( V1 At + uJ2^At - fix) At + o(At)) + 

ox ov V / 

9a d 2 p(x,v,t) .. . 
£1U At dv* +°( At )' 



7 



so 



® gi 



1VCS 



p(x, v , t + At) 



1 -7At 



1 



1 - 7 At 



ax 



At dp(x,v,t) 
1 - 7At 0u 



(«7 - /(^)) 



£7 At d 2 p(x,v,t) . o /9N 



7At 

Dividing by At and taking the limit At — > 0, we obtain the Fokker-Planck equation in the 
form 

d 2 p(x, v, t) 



d ji^A = _ t)M ^) + 8 [(7B _ +£T 
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(15) 



which is the conservation law (HJ) with the flux components 0. The UF Jlr(xi, t) is usually 
defined as the integral of J x (xi,v,t) over the positive velocities 3, and references therein], 
that is, 



JLR{xi,t) 



vp(xx,v, t) dv. 



(16) 



To show that this integral actually represents the probability of the trajectories that move 
from left to right across x% per unit time, we evaluate below the probability flux from a path 
integral. 



III. THE UNIDIRECTIONAL FLUX OF THE LANGEVIN EQUATION 

The instantaneous unidirectional probability flux from left to right at a point X\ is defined 
as the probability per unit time (At) , of Langevin trajectories that are to the left of x\ at time 
t (with any velocity) and propagate to the right of x\ at time t + At (with any velocity), in 
the limit At — > 0. This can be expressed in terms of a path integral on Langevin trajectories 
on the real line as 

^ r x l /'OO POO POO J 

JLn(xi,t) = lim — / d£ I dx drj dv p(£,r},t)6(x - g - rjAt) x 



-777 + /(0] At] 
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Integrating with respect to v eliminates the exponential factor and integration with respect 
to £ fixes £ at x — r]At, so 

Jl p(xi,t) = lim -— / / p(x — nAt,ri,t) dridx 



2 /*0O 

lim — — / dri / 77, t) dw 



r)p(xi,r),t) drj. (18) 
The expression f|T8|) is identical to (fTEj) . 



IV. THE SMOLUCHOWSKI APPROXIMATION TO THE UNIDIRECTIONAL 
CURRENT 

The following calculations were done in ^| and are shown here for completeness. In the 
overdamped regime, as 7 —>■ 00, the Smoluchowski approximation to p(x,v,t) is given by 



. e"^ r , . 1 
p(x,v,t) j= <p(x,t) - - 



dp(x,t) 1 - ( w ' 

— f(x)p{x,t) 

ox e 



v + O 



r 



where the marginal density p(x, t) satisfies the Fokker-Planck-Smoluchowski equation 

dp(x, t) d 2 p(x, t) d 



^ dt " dx 2 
According to flE} and (fl"9"J) . the UF is 



dx 



[f(x)p(x,t)\ . 



JLR{xi,t) 



vp(xi, v, t) dv 



(19) 



(20) 



00 ^—v t2 /2e 



y/2 



p(x,t) 

7T6 I 7 



dp(x, t) 1 
dx e 



f(x)p(x,t) 



v + O 



r 



dv 




e^-mp^t) 



+ o[- 



Similarly, the UF from right to left is 

JKL&ut) = ~ vp(x 1 ,v,t)dv 



(21) 



£ ~dx f( x )p( x ^) 



0[ - 2 



(22) 
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Both UFs in (|21|) and ()22|) are finite and proportional to the marginal density at x%. The 
net flux is the difference 

1 



■Jnet( x i>*) = JLn{x\,t) - JnL(xi,t) 



as in classical diffusion theory Q], [22]. 
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e ?^lL _ f( x )p(x,t) 



(23) 



V. THE UNIDIRECTIONAL CURRENT IN THE SMOLUCHOWSKI EQUATION 



Classical diffusion theory, however, gives a different result. In the overdamped regime 
the Langevin equation (J§J) is reduced to the Smoluchowski equation [^] 

^X = f(x) + y^fW. (24) 

As in Section IIII1 the unidirectional probability current (flux) density at a point x\ can be 



expressed in terms of a path integral as 



JLR(xi,t) = lim J L R(xi,t, At), 



where 



J LR (xi,t,At) 



7 



OO /"OO 



(25) 



(26) 



p (xi, t) — v At 
It was shown in that 

J L R(xi,t, At) = 

Similarly, 



2s 



dx 



°(t) 



ivy At 



27 



<9x 



7' 



jRL(xi,t) = hm^ Jrl(xi, t, At), 



where 



J RL (x 1 ,t,At) 



7 

AneAt ./,, 



°° /"°° f 7<C 2 

d£ J d( exp I — > • 



p(xi,t) + v 7 ^ 



2s 



<9x 



7 



(28) 
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If p(xi,t) > 0, then both JLR,(xi,t) and JRi,(xi,t) are infinite, in contradiction to the 
results (J2U an d (E2J)- However, the net flux density is finite and is given by 



J net( x ^^) = £m {J LR (x!,t, At) - J RL (xi,t, At)} 



1 

7 



' d 

e-^p(xi,t) - f(xi)p(x u t) 



(29) 



which is identical to (|23|). 

The apparent paradox is due to the idealized properties of the Brownian motion. More 
specifically, the trajectories of the Brownian motion, and therefore also the trajectories of 
the solution of eq. (|24|) . are nowhere different iable, so that any trajectory of the Brownian 
motion crosses and recrosses the point x\ infinitely many times in any time interval [t, t + At] 



23| . Obviously, such a vacillation creates infinite UFs. 

Not so for the trajectories of the Langevin equation ((0]). They have finite continuous 
velocities, so that the number of crossing and recrossing is finite. We note that setting 
7 At = 2 in equations flH|) and (J2SJ) gives (jUj) and (1221. 

VI. BROWNIAN SIMULATIONS 

Here we design and analyze a BD simulation of particles diffusing between fixed concen- 
trations. Thus, we consider the free Brownian motion (i.e., / = in eq. (|4*j)) in the interval 
[0,1]. The trajectories were produced as follows: a) According to the dynamics (J3J), new 

trajectories that are started at x(—At) = move to x(0) = x — \Aw\; b) The dynamics 

V 7 
[2e 

progresses according to the Euler scheme x(t + At) = x(t) + * — Aw; c) The trajectory is 

V 7 

terminated if x(t) > 1 or x(t) < 0. The parameters are e — 1, j — 1000, At = 1. We ran 
10,000 random trajectories and constructed the concentration profile by dividing the interval 
into equal parts and recording the time each trajectory spent in each bin prior to termina- 
tion. The results are shown in Figure ^ The simulated concentration profile is linear, but 
for a small depletion layer near the left boundary x = 0, where new particles are injected. 
This is inconsistent with the steady state DE, which predicts a linear concentration profile 
in the entire interval [0, 1]. The discrepancy stems from part (a) of the numerical scheme, 
which assumes that particles enter the simulation interval exactly at x = 0. However, x = 
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/(.<■) C I ex P {-^-^\dy, (30) 



is just an imaginary interface. Had the simulation been run on the entire line, particles 
would hop into the simulation across the imaginary boundary at x = from the l eft, rather 
than exactly at the boundary. This situation is familiar from renewal theory 2M- The 
probability distribution of the distance an entering particle covers, not given its previous 
location, is not normal, but rather it is the residual of the normal distribution, given by 

2a 2 

where a 2 = and C is determined by the normalization condition / f(x)dx = 1. 

7 Jo 

This gives 

«*>=^«*te)- (31) 

Rerunning the simulation with the entrance pdf f(x), we obtained the expected linear 
concentration profile, without any depletion layers (see Figure EJ. Injecting particles ex- 
actly at the boundary makes their first leap into the simulation too large, thus effectively 
decreasing the concentration profile near the boundary. 

Next, we changed the time step At of the simulation, keeping the injection rate of new 
particles constant. The population of trajectories inside the interval was depleted when the 
time step was refined (see Figure EJ). A well behaved numerical simulation is expected to 
converge as the time step is refined, rather than to result in different profiles. This short- 
coming of refining the time step is remedied by replacing the constant rate sources with 
time-step-dependent sources, as predicted by eqs. (|27p -(|28 )) . Figure |U describes the concen- 
tration profiles for three different values of At and source strengths that are proportional to 
l/\/At. The concentration profiles now converge when At — > 0. The key to this remedy is 
the calculation of the UF in diffusion. 



VII. SUMMARY AND DISCUSSION 

Both Einstein ^| and Smoluchowski (see also Q]) pointed out that BD is a valid 
description of diffusion only at times that are not too short. More specifically, the Brow- 
nian approximation to the Langevin equation breaks down at times shorter than I/7, the 
relaxation time of the medium in which the particles diffuse. 

In a BD simulation of a channel the dynamics in the channel region may be much more 
complicated than the dynamics near the interface, somewhere inside the continuum bath, 
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sufficiently far from the channel. Thus the net flux is unknown, while the boundary con- 
centration is known. It follow that the simulation should be run with source strengths (|27jh 

However, Jw is unknown, so neglecting it relative to , / — — r— Cl r will lead to steady state 
boundary concentrations that are close, but not necessarily equal to Cl and Cr. Thus a 
shooting procedure has to be adopted to adjust the boundary fluxes so that the output 
concentrations agree with Cl and Cr, and then the net flux can be readily found. 

According to (|27j) and ([28)1 . the efflux and influx remain finite at the boundaries, and 

' 2 

agree with the fluxes of LD, if the time step in the BD simulation is chosen to be At = — 

7 

near the boundary. If the time step is chosen differently, the fluxes remain finite, but the 
simulation does not recover the UF of LD. At points away from the boundary, where correct 
UFs do not have to be recovered, the simulation can proceed in coarser time steps. 

The above analysis can be generalized to higher dimensions. In three dimensions the 
normal component of the UF vector at a point cc on a given smooth surface represents 
the number of trajectories that cross the surface from one side to the other, per unit area 
at x in unit time. Particles cross the interface in one direction if their velocity satisfies 
v ■ n(x) > 0, where n(x) is the unit normal vector to the surface at x, thus defining the 
domain of integration for eq.(|5|). 

The time course of injection of particles into a BD simulation can be chosen with any 
inter-injection probability density, as long as the mean time between injections is chosen so 
that the source strength is as indicated in (|2*7j) and (|28)1 . For example, these times can be 
chosen independently of each other, without creating spurious boundary layers. 
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Concentration vs. Distance: Equal Source Strengths 
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A. Singer, PRE, Fig. 3 
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Concentration vs. Distance: Equal Source Strengths 
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FIG. 3: The concentration profile of Brownian trajectories that are initiated at x = and termi- 
nated at either x = or x = 1. Three different time steps (At = 4, 1,0.25) were used, but the 
injection rate of new particles remained constant. Refining the time step decreases the concentra- 
tion profile. 
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Concentration vs. Distance: 1/sqrt(dt) Source Strengths 
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Concentration vs. Distance: 1/sqrt(dt) Source Strengths 
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FIG. 4: The concentration profile of Brownian trajectories that are initiated at x = and termi- 
nated at either x = or x = 1. Three different time steps (At = 4,1,0.25) are shown, and the 
injection rate of new particles is proportional to 1/y/At. Refining the time step does not alter the 
concentration profile. 
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